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ABSTRACT 

Recent interest in military VTOL/STOL aircraft employing unprepared 
landing sites has led to interest in the problem of landing surface ero- 
Sion. Surface erosion is caused by the aerodynamic forces on ground 
particles existing within the flow field of an impinging jet. The 
inviscid flow field is discussed and the viscous ground boundary layer 
is analyzed utilizing both theory and available experimental data. A 
mathematical model of the process of entrainment of ground particles is 
constructed. Erosion rates in the form of erosion profiles are pre- 
dicted for selected jet configuration and types of terrain. A criterion 
for entrainment, due to both lift and drag, was found and presented for 


selected distances from the jet centerline. 
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NOTATION 
erosion rate - ft/sec 
coefficient of drag 
coefficient of skin friction 
nozzle diameter - ft. 
particle diameter - in. 
force - Lb. 
lift force due to the effect of the ground - lb. 
lift force due to shear - lb. 


mass - lb sec“/f£t. 


rotational flow parameter = ae r 
ele 


Reynolds number = LS 
I 


number of strips used in strip analysis 

Static pressure on the ground lb/ft. 

shear velocity - ft/sec. 

radial distance from jet centerline - ft. 
radius of particle - in. 

kinetic energy - ft Lb. 

velocity parallel to the ground - ft/sec. 
velocity perpendicular to the ground - ft/sec. 
coordinate parallel to the ground 

coordinate perpendicular to the ground 


vertical distance at which VU = Umax in turbulent boundary layer 
ez 


height of nozzle above the ground - ft. 


density - lb sec2/£t*, 


viscosity - lb sec/ft*, 





=, angle measured from center of particle 
Ho vortic Lty 


5 boundary layer thickness 


Subscripts 


crit. value at boundary layer transition 


J indicates value associated with jet 

L value at local point on the particle 

a property of air 

I value in inviscid flow 

P value associated with a particle 

Cc measurements made on cylinder in strip analysis 
S static 

A ambient 

n measurement made on a strip 

W value associated with effect of wall constraint 
S value associated with shear effect 


vi 





l. Introduction. 

The increasing interest in, and use of, VTOL aircraft and other 
vehicles with vertically directed fans and jets, particularly in mili- 
tary operations, has brought about concern for the deterioration of 
unprepared or soft landing sites. The downwash impingement caused by 
these aircraft results in a radially developed viscous flow field on the 
ground, producing forces adequate to cause entrainment of the ground 
particles. The reasons for concern with this problem are numerous. [In 
addition to possible damage to the aircraft itself and the erosion of 
the landing surface there are the accompanying hazards to ground per- 
sonnel, damage to ground support equipment, poor visibility afforded the 
pilot, and the violation of the military concept of concealment. The 
present study was designed to construct a mathematical model of the 
ground erosion process to allow prediction of the initial erosion rates 
for various types of terrain for various nozzle velocities and geome- 
tries. Previous work ‘Ee iG. 154. 22,5. 23, 25 | on the problem of erosion 
due to downwash impingement has been primarily an experimental examina- 
tion of the problem with small scale jets. 

The critical mechanism in the erosion process is the interaction of 
the ground boundary layer and the terrain particles. An analysis of 
this mechanism requires an understanding of the viscous mixing problem 
of a finite free jet, with the inviscid flow field extending radially 
from the jet center line as a result of the presence of the ground plane. 
An analysis of the phenomenon which produces entrainment of the ground 
particles is basically a study of the aerodynamic forces which exist in 
the ground boundary layer, as they are affected both by the decaying in- 


viscid flow field above the boundary layer and by the constraint of the 





ground surface. In analyzing the boundary layer forces available theory 
was utilized whenever possible, but since the radially developed bound- 
ary layer runs the gamut of stagnation point flow to turbulent flow and 
then complete decay, there are areas for which no theory or exact solu- 
tion exists. To fill these holes, experimental data were used. Where 
exact solutions were used for the velocity profiles within the boundary 
layer, verification was made with experimental results if data existed. 
The analysis was restricted to single jets of uniform dynamic pres- 
sure loading and circular cross-section under "no wind'' conditions shown 
diagrammatically in Fig. 1. An axi-symmetric boundary layer was assumed, 
acting over a terrain of uniform particle size and density. Table [ in- 
dicates the various combinations of nozzle characteristics and radial 
stations that were analyzed. The prediction of erosion rates is neces- 
sarily restricted to incipient motion, since the secondary effect of 
random collision is not readily treated by analysis and because the con- 
stantly changing shape of the ground surface is not accounted for as a 
continuing process. Although this work is restricted to consideration 
of uniform jets only, previous experiments by Morse (23] indicate that 
the dynamic pressure distribution due to downwash from rotor blades and 
ducted fans is similar to that for a uniform jet. The uniformity of 
the jet or downwash appears to be a critical factor only for ratios of 


jet altitude to diameter, Z/Dy, less than unity. 





Zi History of the Problem. 

The classic problem of an impinging jet on a normal surface is not 
new, but its application to the field of VTOL aircraft, of either rotor 
or jet type, has been generated only within the past several years. The 
problem of surface erosion resulting from impinging downwash first re- 
ceived attention from Kuhn hi4 Kuhn's work was experimental employing 
small scale test equipment to make dynamic pressure distribution surveys. 
An experiment was conducted by fixing the nozzle dynamic pressure and 
raising or lowering it over terrain of known condition. Testing was 
done with various types of sod, dry and wet dirt, dust, and sand. Simi- 
lar experimental tests have been conducted more recently by Morse \.231\ 

The most recent work done on this problem is that of Brady and Lud- 
wig (26. This work consists of extensive experiments to investigate 
the characteristics of the impinging uniform jet, including both the 
inviscid flow and the ground boundary layer. Measurements of boundary 
layer velocity profiles were compared with theory and showed agreement 
to a limited degree. Some of these experimental results have been used 
in support of the present work. Vidal (28 | has given an insight into 


the aerodynamic forces existing within the boundary layer. 





or The Flow Field. 

In order to consider the mechanism of particle entrainment and ground 
erosion, a thorough knowledge of the flow field of the impinging jet is 
required. This flow field can effectively be separated into three re- 
gions. The first of these regions exists regardless of the presence of 
the ground surface, and is the region of viscous decay immediately exter- 
nal to the jet nozzle. Viscous decay of a free jet is a classical prob- 
lem under uniform conditions. This region is characterized by the flar- 
ing out of the jet as the jet boundary mixes viscously with the stationary 
field external to it. Introducing a ground surface into the flow field, 
normal to the jet, does not invalidate the quantitative solution of free 
jet decay, but restricts its extent of usefulness. Experimentation has 
shown that the theory of viscous decay of a free jet is adequate for 
the impinging jet problem for nozzle heights of greater than about eight 
jet diameters. Since in this work, the interest is in nozzle heights 
of Z/f,, << g » an additional region of viscous mixing with the station- 
ary boundary must be studied and determined, beginning with the point at 
which the presence of the ground surface has altered the free jet charac- 
teristics. Further study must then be made of the viscous decay of the 
jet after impingement on the ground plane. To date no satisfactory 
theoretical method exists for dealing with the entire flow region in 
three dimensions. This problem may be partially overcome, however, by 
including the viscous decay region in describing the inviscid flow field. 

Knowledge of the inviscid flow field is necessitated by the require- 
ment that the inviscid velocity existing at the upper edge of the ground 
boundary layer be known. The difficulty in finding a theoretical solu- 


tion for this region lies in the fact that although the boundary conditions 


are well defined, the location of the free streamline is not known. Three 
dimensional theoretical solutions of this problem are limited to approxi- 
mate methods. One of these, by LeClerc {10], uses an electrolytic analog 
to determine the shape of the free streamline boundary. In this method 
the solution of the inviscid flow region was accomplished by relaxation 
techniques, predicting velocity and pressure distribution. Such a 

method of solution could be made to approach the exact solution to any 
desired degree of accuracy. 

In the present work an exact solution of the inviscid flow field 
was not required. [It was assumed adequate to use existing experimental 
data which compared favorably with the limited theoretical analyses avail- 
able. In most instances it was found that the theoretical solutions were 
not in good agreement with the experimental data for oy <4 - The approx- 
imate methods of predicting the characteristics of the tavieeta flow re- 
gion by use of an equivalent inviscid jet of greater diameter and decreas- 
ing dynamic pressure, or by replacing the jet with a cylindrical vortex 
sheet of constant radius extending to the ground, have proven unsatis- 
factory when compared with experimental data. In view of this, available 
experimental data were used in this work as a solution to the inviscid 
flow field. Fig. 2 shows the result of empirical equations fitted to 
these data. 

The third flow region is that of the ground boundary layer. A thor- 
ough knowledge of this region was most important to this work as the ground 
particles are predominately immersed in the boundary layer, and the mechan- 
ics of particle movement originate with the aerodynamic forces resulting 


from the characteristics of the boundary layer. 


4, The Ground Boundary Layer. 

The key to the entire analysis of ground erosion is the mathematical 
model used for the boundary layer. Sound theoretical analysis of the 
boundary layer was used in this study when possible. Experimental 
results were used whenever the theory was non-existent or was in large 
disagreement with these results. For effective analysis the boundary 
layer may be divided into four regions. The first is the stagnation 
area boundary layer existing adjacent to the jet center line and extend- 
ing to Rip, = -2 or to a station directly under the original free jet 
boundary. This area is laminar for the conditions investigated. The 
solution of this region utilized the Himenez solution (18 | for the uniform 
impinging jet on a perpendicular wall. Experimental results were not 
available for this region, so that the validity of this solution is not 
verified but was assumed. The Himenez solution in three dimensional 
flow, developed by Homann [18 | is an exact solution utilizing the Navier- 
Stokes equations for axisymmetric flows. The solution was developed 
explicityly for stagnation conditions and the velocity profiles result- 
ing from this have been employed to Rip y= 0.5. The nondimensional 
velocity profile is shown in Fig. 3. 

The second region is the laminar boundary layer extending from the 
edge of the original jet boundary to the radial station at which transi- 
tion begins. A complete method of solution for this laminar region, by 
Smith [19], is well verified by experiment for +/p,,= .5 \26|, but 
gives little insight as to the limit of the laminar region and the be- 
ginning of transition. This point has been approximated by Brady and 
Ludwig [26] from considerations of neutral stability of the laminar bound- 


ary layer. A Reynolds Number at which the boundary layer becomes neutrally 


stable can be determined [18 | and may be converted to a jet nozzle 
Reynolds Number as a function of PJon) crit: On this basis the 
boundary layer becomes unstable at oy SF 1.0 for the jet velocities 
and nozzle diameters analyzed in this work. In view of this, and the 
fact that the invisced velocity generally reaches a peak in the vicinity 
of Riso = 1.0 (Fig. 2), the beginning of transition was taken to be 
a5, =O. 

The transition region is even more difficult to define. Since 
fully turbulent flow exists theoretically when the pressure gradient on 
the ground surface is essentially zero, it was assumed that the transi- 
tion region extended to the station where the pressure gradient could be 
taken to be negligible. Fig. 4 shows the static pressure distribution 
near the ground taken from an electrolytic 3 of the entire flow 


field [10]. The station at which the gradient, ; EN, Ree) / 8 a 


was taken to be zero was Ryo, = 2.0. For the laminar and transition 
regions of the boundary layer experimental results were used to form the 
velocity profiles. This data, collected by Brady and Ludwig [26 | for a 
limited range of jet velocities indicate that the non-dimensional 
velocity profile was almost completely independent of the jet velocity. 
The measurements were taken at nozzle heights in the range Z/Oqy= 0.5 
to 4.0 at four locations and for the radial range of y/o te to 1.33 
at six different stations. These twenty-four configurations are assumed 
to be fairly representative of the laminar and transition regions and 
are shown in Figs. 5a through 5h. For purposes of computation, curves 
were fitted to each group of data. 


The fourth region in the boundary layer is the turbulent regime. 


It was assumed that this region begins at Sony 2 as this is the 





Station at which the pressure gradient is assumed negligible. A Limited 
amount of work has been done on theoretical solutions of turbulent bound- 
ary layers. One of these that appears applicable to the impinging jet 

is that of Glauert [7]. This single solution for velocity profiles in 
the turbulent boundary layer, assumed valid to Rio w= 4. which is the 
radial limit of the analysis, has been verified by experimental results 
at R/ay = 1.33 [26]. The uncertainty as to the characteristics of flow 
decay limit further extension of the analysis. 

Glauert set up the boundary layer equations and found a similarity 
solution in which the form of the velocity distribution across the jet 
does not vary along its length. The solution is dependent upon the 
assumption of an effective eddy viscosity as required to satisfy the law 
due to Blasius for flow in a pipe as well as Prandtl's hypothesis for 
free turbulent flow. The first solution is considered to be valid near 
the ground surface, and the second to be valid above some determined dis- 
tance from the ground surface. Fig. 6 shows these results in the form 
of a velocity profile, for the radial distances and nozzle heights con- 
sidered in the present work. It was found that the shape of the velocity 
profile was essentially independent of both the radial distance and the 
nozzle height, but is dependent upon the nozzle Reynolds Number. 

The turbulent region is characterized by a velocity profile which 
approaches a maximum, U moy , with increasing y and then decreases with 
further increase in y. As is shown in Fig. 1 it is expected that the 
turbulent boundary layer will continue to grow as long as there exists 
an inviscid upper boundary. Since very little attention has been given 
to this phenomenon of flow decay, and because of the nature of the tur- 


bulent solution, a fictitious boundary layer thickness is assumed for 





this region. This is defined to be the thickness that exists at 
U=Umoy + It was found that the nature of this assumption does not 


affect the overall results appreciably. 





ee Analysis of Aerodynamic Forces Within the Boundary Layer. 

The ground boundary layer is a non-uniform flow that is character- 
ized by a large velocity gradient or shear layer extending over a large 
portion of the thickness. In addition, the ground surface or wall bound- 
ary at the base of the flow forms a constraint on the streamlines as 
they attempt to pass beneath an obstacle bedded in the ground surface. 

There are three distinct forces that exist within the boundary 
layer. These are: 1) Drag due to the horizontal velocity component of 
the flow; 2) Lift due to the proximity of the wall; and 3) Lift due to 
the velocity gradient in the flow. There is no theory available which 
deals with this problem in general. Therefore it was necessary to piece 
together an approximate theory for each of the various forces. 

Drag Force. 

The drag force was the easiest to account for. To simplify the 
analysis in a non-uniform flow, a strip analysis was made in which the 
spherical particle was divided into a number of layers and it was assumed 
that a uniform flow acted upon each layer. Knowing the velocity of the 
flow on each layer, from the velocity profile, the force on each layer 
and the summation of horizontal forces were developed. The calculation 
of the total drag force employed the drag coefficient of cylindrical 


bodies measured by Hoerner |9 |: 


Nir a. YOO 


ec or 
Cy = i 


100 < Ne < 10% 
=. ae 
an = Loe 


10 





In addition to the drag force the skin friction force was calcu- 
lated for each layer and summed for the particle. The friction force, a 
result of a minute boundary layer effect on the particle itself, was 
relatively small but in some cases significant. For spheres in low Rey- 


nolds Number flows the skin friction coefficient can be approximated by: 


Cre = C6 FN R)/2 


The tendency of the velocity gradient to cause a moment on the particle 
was neglected. 
Lift Due to Wall Proximity. 

Theories accounting for the effect of a wall constraint in non- 
uniform flow do not exist, but this effect was estimated using the exist- 
ing theory for uniform flow [12]. The effect of the ground surface con- 
straint is to distort the flow over the sphere, crowding the streamlines, 
causing a negative lift force. The theory is based on kinetic energy 
considerations, maintaining equilibrium on a sphere in the presence of a 
wall. For a sphere, in a flow parallel to the wall, an upward force is 
required to maintain equilibrium, indicating that the sphere must be 
attracted to the wall. 

The method of calculation of this lift force was similar to that of 
drag. The strip approximation method was used to determine an average 


local velocity over a particular strip. It was then assumed that this 


average local velocity was that existing in the vicinity of the stagnation 


ll 





point, realizing that the stagnation point is shifted toward the region 
of higher velocities. The lift then becomes a function of the sphere 
size and the average local velocity, for a particle bedded on the ground 
plane. The kinetic energy of a sphere in a moving fluid near an in- 
finite fixed wall can be approximated by a first order solution [12}: 


T=’ CAV + BU” ) 


where: 


D 
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Mo. = mass of fluid displaced by the sphere. 


From Lagrange's equation for equilibrium; 
a ee eee a = 
de Cx) - rx =O 


If the sphere is to be maintained in equilibrium, the result is a force 


exerted on the sphere; 
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Making the substitutions, the force exerted on the sphere in the y or 


vertical direction becomes: 
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where: 
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Mo = <6. 730F 

we get; 
a a 6 

a = 6, < UD, ‘na c A+ 

Assuming that the particle is very close to the ground such that; 


y-r<<l 
or 


the Lift due to the wall constraint becomes: 


Ly ~ Veta geo ae 


For the strip analysis, summation yields: 


a a 
lw = Fie Ra ny 2 On) fp 


Lift Due to Velocity Gradient. 

There have been a number of theoretical studies of the effects of 
velocity gradients or shear on aerodynamic forces. Most of these, how- 
ever, deal either with an unbounded flow or two dimensional flow. [In 


comparing two dimensional solutions with three dimensional solutions it 
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was noted that the difference is a term in the three dimensional solu- 
tion describing the lateral component of flow. For the axi-symmetrical 
case dealt with here, this component was assumed negligible. A solution 
for non-uniform flow on cylinders by Murray and Mitchell [13] was em- 
ployed in conjunction with a strip analysis. 

In the strip analysis, the particle was divided into layers of cir- 
cular cylinders situated such that their centers formed a line perpendicu- 
lar to the flow. Beginning with the final result of Murray and Mitchell, 
the dimensionless shear velocity on each cylinder was found to be: 
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YCID = Euler's Constant 
= 0.57721567 


The dimensional shear velocity is given by: 


$= 9 Ue 
which is of the form (when squared): 
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(J, is the local velocity at the center of each cylinder, and I, K and 
K' are Bessel functions of the first and second kind, and the derivative 
of the second kind respectively. N is a flow parameter dependent upon 
local velocity, radius of the cylinder and the vorticity at the stagna- 
tion point, and is defined as: N= 6 Ke . The determination of 
the flow parameter was approximated since the stagnation point on each 
cylinder could not be known in advance of the calculations. In order to 
keep the shear velocity in the direction it was known to flow, a lower 
bound of N=] could be established. If “Jo and U are evaluated at 
the geometric centers of the cylinders, experimentation has indicated an 
upper bound of N= 3 ia3le After analyzing the computed results of 
shear velocity with values of N ranging from 1.1 to 10, a value of 
N = 1.5 was accepted as being representative of the physical model. 
Having found the shear velocity the lifting force due to shear 
effect was calculated for each cylinder and summed over the particle. 


The resulting expression is: 
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Where: 
radius of cylindrical strips 


i = 
b = width of strip 
Yr. = number of strips 


Noting that: 
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In evaluating the integral we have: 
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Further we may write: 


Since the last two integrals are identically zero. 
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Since the first and last integrals are identically zero. 

Due to the complexity of the above expressions and the difficulty 
of integration, it was desired to check this method with another method 
29] of a much more approximate nature for calculating the lift due to 
shear. 

The alternate method of calculating lift assumes that the previously 
discussed lift due to the ground surface constraint is negligible and 
that the lift forces on the particle can be estimated by using the solu- 
tion of Hall 8 for a sphere in a two dimensional uniform shear flow. 

In addition, the approximate method assumes that the sphere rests on a 
bed of spheres of the same diameter and that the pressure on the bed is 
the ambient stream static pressure. The lift coefficient then becomes 
simply a function of the local slope of the velocity profile, ; 
the local velocity on the sphere, and the radius of the sphere. In view 
of the many restrictions placed on this solution it was used only as an 
order of magnitude check on the other solutions. Fig. 7 shows a compari- 


son of the two solutions. 
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6. Particle Entrainment and Incipient Erosion Rates. 

Only a few of the possible types of terrain to which a mathematical 
analysis of the erosion rates might be practically applied, were investi- 
gated. Uniform particles ranging in size from .003 in. to .125 in. in 
diameter were investigated. The smallest sizes represent loose dirt 
and the larger diameters represent sand or sandy gravel, with .125 in. 
being representative of small gravel type rock. It was assumed that all 
terrain was without any moisture content. The densities of the particles 
ranged from 75 lb/£t3 for loose dirt to 125 1b/ft? for gravel. Table II 
lists the various terrain particles analyzed. 

There are three possibilities for entrainment of the ground particle. 
The first is that it may have enough lift force on it to be lifted directly 
from the ground surface into the free stream. The second is that after 
rolling motion is started due to the drag force, it may be able to gain 
the additional lift required for entrainment in the free stream. The 
third is that with many particles rolling at different velocities the 
random collision of particles that follows will bounce some of them into 
the inviscid stream. Only incipient motion is within the scope of this 
work, hence it was assumed that any particle which was moved contributed 
to the erosion rate. 

The net drag force and the net lift force, including both mechanisms 
of lift, were computed for each of the particles, at each of thirteen 
radial stations on the ground surface for each of the nozzle configura- 
tions studied. To calculate the initial net drag force, the previously 
mentioned ground surface model of nested particles was assumed using a 


coefficient of static friction of ; 
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Le Results and Conclusions. 

The particle size criteria for lifting entrainment and drag entrain- 
ment are shown in Figs. 8, 9, and 10 in which the particle size is plotted 
versus a load parameter at various stations with jet velocity and jet 
height as parameters. The load parameter is defined as the particle den- 
sity divided by the dynamic pressure in the free stream. The Lift en- 
trainment criterion predicts that an optimum particle size, the particle 
size most prone to entrainment, exists for all radial stations and that 
the critical density increases with radial distance to a certain Se 
The optimum particle size for lift entrainment occurs because particles 
smaller than the optimum are affected by much smaller velocities, and 
particles larger than optimum are affected more by a flow region where 
the shear gradient, mo , approaches zero. 

The drag criterion indicates that extremely small particles of 0.02 
in. diameter or less, up to very high densities will roll on the ground 
(see Fig. 8). As the particle size is increased the possibility of 
rolling decreases rapidly up to a particle diameter of 0.03 in. at which 
size the tendency to roll increases again, since the larger velocities 
in the boundary layer begin to act on the particle, until an optimum 
diameter of approximately 0.06 in. is reached. For radial distances 
less than Rlon = ~l1, velocities in the boundary layer are not suffi- 
cient to move large particles of high density. The range of optimum size 
particle for drag entrainment of approximately .03 to .09 is consistent 
for all radial stations. It is of particular significance that the range 
of optimum particle diameter for both lift and drag is very nearly equal 
to the boundary layer thickness at the particular radial station being 


considered. 
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The incipient erosion rates for selected densities, particle sizes, 
and nozzle configurations are shown in Figs. 1l and 12. These are found 
by combining as vectors the incipient horizontal and vertical accelera- 
tions of the particles. This was assumed to be proportional to the 
incipient erosion rate. The erosion rate, normalized by the maximum 
rate, is plotted versus the nondimensional radial distance. For all con- 
figurations and particles the erosion profile is maximum in the area 
oa ye <2.0. The maximum erosion for larger particles occurs near 
Voy = 1, and then drops off sharply with increased radial distance. 
As the particles become smaller a maximum erosion occurs closer to the 
perimeter of the jet nozzle and a second maximum begins to occur at 
Ring = 1.3. This occurs after the transition region has begun to 
form and where the onset of turbulence can affect the particles. Fig. 13 
shows the incipient rate of lift entrainment by itself and indicates that, 
in general, the shape of the erosion profile can be attributed to the 
lift forces, 

The results shown in the graphs represent only a few configurations 
selected from the large amount of data collected. In viewing all of the 
results and comparing the many combinations of configurations, it was 
found that the variation of ground erosion rates was principally a func- 
tion of the dimensionless radial distance from the jet centerline and 
the size of the particle for a given nozzle height, ye . The ero- 
Sion rate was nearly proportional to the particle density except near the 
stagnation point, where the erosion rates were very small. 

From these results it may be concluded that for ground particles of 


less than .04 in. diameter the incipient erosion rate has two maximum 


points. One in the vicinity of the perimeter of the jet nozzle, and the 
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second at a radial distance approximately two diameters from the stagna- 
tion point. For very small particles and low jet velocities the second 
maximum is the most significant. [In ali other cases where two maximum 
points exist, the first is the most significant. For ground particles 
whose diameters are greater than .04 in. a single maximum point exists 
near the perimeter of the jet nozzle. As the particle size increases 
the influence of the turbulent region decreases. 
It was concluded that particles of approximately .02 in. diameter 
or less would entrain due to drag even at very high densities or very 
low velocities. For particles larger than approximately .02 in., an 
optimum size of approximately .06 in. exists up to radial distances of 
four nozzle diameters for drag entrainment. An optimum particle size 
most susceptible to entrainment due to lift exists in the range of .0/ 
to .12 in. diameter for radial distances up to four nozzle diameters. 
Also, it may be concluded that varying the nozzle diameter had very 
little effect on the dimensionless erosion rate; hence in the final analy- 
sis the nozzle diameter was not considered as a parameter. Moreover, 
while the actual erosion rates were approximately linearly dependent upon 
the jet velocity, the dimensionless erosion profile appeared to be inde- 
pendent of jet velocity, with the exception of the order of importance 


of the maximum points. 
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APPENDIX 


FORTRAN PROGRAMS 


LIFT FORCES: 


Input. 


Program is generated from the polynomials describing 


the velocity profiles of the boundary layer, for 


various radial stations. 


MM 
PCOZ 
QSH 
Vd 
DN 
RG 
2G 
YHVX 


vMXx 


Number of terms of the polynomial 

Coefficients of the polynomial 

Nondimensional shear velocity 

Jet velocity 

Nozzle diameter 

Radial distance from jet centerline 

Nozzle height above ground 

Distance to have maximum velocity in turbulent 
boundary layer 

Maximum velocity in turbulent boundary layer 


intermediate Results.and Output 


a 
Vi 
YND 
XND 


DNS 

DVX 

SDVX 
DUDY 
F'LSG 
SLFT 
PLET 
PLAC 


Particle diameter 

inviscid velocity 

Nondimensional distance above ground 
Nondimensional roots of velocity profile 
polynomial 

Particle density 

Derivative of polynomial 

Slope of velocity profile - 

Slope of velocity profike - 

Shear lift (cylinder method) 

Shear lift (sphere method) 

Lift due to wall constraint 

Acceleration of particlo in vertical direction 


Subroutine Root finds the roots of tne velocity profile 
polynomial 


pe 





SUBROUTINE LIFT( MM,PCOZ» QSH,»,VJ»DN,RG»ZG,NCC, YHVXs VMX) 
DIMENSION DP(10),RP(10), YGU(40),YGB(40),YL(70).YG(70), 
4AR(79)+VL(70)-PCO2(350),XL(490),FX( S50) ,NLY(5),DVX(50), 
2 ce erectile aa os 
S Vili) sy oDuU Ceo 
DP(1)=.06 
pete) =.04 
He os 02 
DRt4)=.01 
hie Oren 7 
DP(6)=.125 
pP(7)=.0038 
RGN=RG/DN 
ZGN=ZG7DN 
PRINT 641, DN 
6141 FORMAT(/S// 2ZHNOZZLE DIAMETER (FT) = F9-2 ) 
PRINT 612, VJ 
612 FORMAT( 26HNOZZLE VELOCITY (FT/SEC) = F7.4 ) 
PRINT 400,RGN 
400 FORMAT(C26HNONDIMEN RADIAL DISTANCE = F20.10) 
PRINT 401, ZGN 
401 FORMAT(28HNONDIMEN HEIGHTH OF NOZZLE = F20.10) 
DO 245) 1isis7 
PRINT 40%, DP(11) 
403 FORMAT(//25HPARTICLE DIAMETER(IN.)= Fe0.10): 
RP(IIT)=DPCII)Z2. 
GNU=S3.745E-7/2.3769E-3 
IFCNCC-25) 502, $03, 592 
303 vil=VMX 
GO To 700 
302 IF C(ZG/DN-.5)640,640,661 
661 KZG=ZG/UN 
GO 10(640,641,640,643),K2ZG 
640 IF CRGN-.$5)60,60,61 
61 IF (RGN-.75)18.18,19 
19 IF(RGN-1.2)20,20,22 
60 VI=ZVU*¥SQRTF CABSF(1.88*RGNe*2-,225*RGN) ) 
GO. 10 od 
18 V]=VU*SOKTF(-3, Z2aRGNwee+d, 15 eRGN-1.26) 
Gor 10 Via 
20 VI=ZVIU*SQRTF(-2,68eRGNee2+9, 48 *RGN- 1.82) 
60 TO 706 
22 VI=VU*SQRTF(2.25*RGNe* 2-9, 45eRGN+8B,96) 
GO TO 700 
644 IF(RGN-.35)6q,60,644 
644 IF (RGN-1.0)645,645,647 
645 VI=ZVI*¥SORTFE (-.0985*RGN¥*¥241.135*RGN-.241) 
GO 10 700 
647 VISVUeSORTFE(.077e#RGNew2- 6491 eRGN41, 599) 
GO To 700 
643 IF (RGN-.40)62,62,654 
654 1F (RGN-1.2)655,655,65/7 
655 V1]=VU*SORTF(-.5eRGNe weet. 2eRGN-.S) 


92 








GOT Ong 7el0 
$57 VI=VU*eSORTR(.0409*#RGNe*2- , S18 *RGN+. 744) 
GO TO 700 
62 VI=VI*SORTE (. 62D we RGNwe wet, Q2>eRGN) 
700 ANCL=101. 
916 TH=DP(II)/ ANCL 
NN=ZANCL 
NA=(NN*4)0/2 
R(1)=RPCII) 
pie 4. K=2,NA 
AA=K-1 
XLC=AA*TH 
4 ROCK) =ZSORIFCRPC ILI) *#*2-XLC we?) 
100 IF (NCC-26)901,900,900 
301 IF (NCC-25)704,705,900 
705 YND( 1)=R(1)/YRVX 
GO: 10 104 
704 YND( 1)=(R¢ 1)/612.*DN))*SORTF(VU*DN/GNU) 
GO TO 104 | 
7090 AAAH=VISRG 
YND( 1)2R(4)*«SQRTF(2.*AAA/GNU)D/12. 
104 XND(1)=0.0 
IF (NCC~25)508,504,508 
504 IFCYND( 1)-.15)500 .500,501 
500 MM=8 
DO 541 I1Z2=1,.MM 
ta PCOX C2) =PCco7C]Z) 


GO TO 50a 
301 MM=4 
DO 202 y=1,™M™ 
N=] +8 ° 


502 PCOX(1)=PCOZ(N) 
506 DO 507 Kw=1.5MM 
EXL=MM-Kw 
FXCKW)I=PCOX(KW)*YND( 1) *wEX4L 
507 XNDC 1)=XNDO 1)4*FXCKW) 
GO TO 673 
508 DO 71 KZ=1.™MM 
EX1=MM -KZ+1 
FX(KZ)=PCOZ( KZ) *YND (1) EX] 
714 XND(1)=XNDO1)+FXOKZ) 
GO TO 67% 
EUS (GALL ROO! (Ll aMNeRCOZS NOS NDSeANIEC ) 
SS ND Te 0) Sis Ola oe 
32 XND(1)=.9999 
31 VL(1)=V1*XND(1) 
706 ART=0. 
DO 43 JZ=1,NA 
Percd2-1 344544545 
44 ARE=2.*R( UZ) *5.149159%TH 
GO TO 43% 
45 ARE=4.¥*R(UZ)43,149159*"TH 
4$ ARTZEART+ARE 
SPA=RPC IIT) *¥*2"4,4%3,.14159 
ADS=SPA~-AR} 
ADS=ADS/SPA 
ACOR=1.+ADS 
reS=6:. 
DO 39 IZ=1,NA 
Pe Clee aes Oy 4d 
40 FLE-TH/2.*2.35769E-S¥R( TZ) eVL(1L)**2eQSH/1494., 
GO 10 39 


ao 





FLSC=FLS*ACOR 
APPROXIMATE SHEAR LIFT OF PARTICLE BY SPHERE METHOD 
V0 Velvet eR 1) 
PROONCC=2509 S2777+129 
77 XNNODD( 1) sYGDUt tT) YR x 
BO Ua oO 
129 YNDD(1)=YGDD(1)*SORTF(2.*AAA/GNU) /12. 
GOO. Gs 
132 YNDD(1)=(YGDN(1)7(12.*DN)) *SQRTF(VU*DN/GNU) 
133 SDVX=0.0 
DO 74 KD=1,M™™M 
MX2=MM-KD 
EXL=MM-KD+1 
fPeONGOC 22) / 257957 5 
79 PCOZ(KD)=PCOX(KD) 
75 JTF CKD-MM)180,181,180 
180 DVX(KD)=EX1*PCOZ( KD) «YNDD(1) **MX2 
Con 10 44 
181 DVX(KD)=PCOZ(KD) 
74 SDVX=SDVX+DVXC(KD) 
IF (SOVX)182,185,188 
162 1F(NCC=25)184,1655184 
184 SDVX=0.0 
IF (NCC -25)190,183,192 
183 DUDY=SDVX*®VMX/YHVX*12. 
GO TO 196 
192 DUDY=SDVX*SOGRTF (2. *AAA/GNU) * VI 
GO TO 193 
190 DUDY=SDVX*SQRTF (VU*DN/GNU) *®VI/DN 
TOS 7 ON = RC LO Vl ey Ze 
CLFH0-345+0-557¥Z0N+0.8/9«%Z0N¥*2 
SLFT=CLF*2.3769E-3S¥VL(1) ee 2/2. "5.149159 HR(1) 2/1494, 
PRINT 800 


800 FORMAT(14HNONDIMEN SHEAR,5X,17HNONDIMEN VELOCITY,5X,17HLIFT DUE TO 
1 SHEAR,5X,17HLIFT DUE TO SHEAR,®X,19HAREA DIFFERENCE 0/0,5X,15HPRO 


ere SeOre, 
PRINT 801 


801 FORMAT(3X,8HVELOCITY,11X,11HON CYLINDER,» 12X,10H(CYLINDER),435X, 


1 8H(SPHERE)/) 
PRINT 602,GSH,XND(1),FLSC,SLFT,ADS,SDVX 
802 FORMAT( 6620.10) 
COs 52 
5 PRINT 200 
200 FORMAT (/5HERROR/) 
52 NNI=59 
UN=NNT 
603 pGBi1)=0. 
YH=0. 
YLM=NNT-1 
YINC=DPCII)/YLM 
DO 169 J=1,NNT 
Pe CdR 7666; 111 ¢ii2 
4111 YG(J)=0. 
GO TO 169 
112 YGC(I0-=YH*YING 
Vina Yd) 
169 CONTINUE 
DO 115 K=2,NNT. 
YND(10=0. 
IF (NCC-26)116,902,902 


OE a 





902 


eS 
$01 


300 
113 
671 


720 
724 
B11 


702 


702 


708 


37 
672 


YNDCKI=YGCRI*SQRTF (2. *AAA/GNU)D 12. 
GO TO 116 
Lene eo) 500s SUT oc 
YNDCK)=YGCKI/SYHVX 

GO To i138 
YNDCKIECYG(K)/(12.*0DN)) *SQRTF(CVJ*DN/SGNU) 
CONTINUE 

DO 672 NX=14,NNT 

XNDCNX)=0. 

JF (NCC=+25)708,720,708 
PEGCUNDUNKI= 379) 721 6 74197 01 
MM=8 

DO 811 14=1,MM 
PCOX(1Z)=PCOZ( JZ) 

GO TO 708 

MM=4 

DO 702 [=1.™MM 

N=1+8 

PCOX(I)=PCOZ(N) 

DO 707 Kw=i5MM 

EXL=MM-KW 
FXCKW)=PCOX( KW) *¥YNDCNX) ®* EX1 
XNDONX)=XNDONX)4+FXCKW) 

GO TO 672 

DO $7 KL=15MM 

EXI=MM-KL+1 
FXCKLISPCOZ( KL) *YNDONX) *®wEX1 
XNDONX)=XNDONX)+FXCKL) 
CONTINUE 

DO 95 KC=1, NWT 
EECAND Cn.) 16 $3:55.5.55 

IF (KC0-1)058,38,54 

XNDCKC)=0.0 

GO 10355 

XNDCKCI=EXNDCKCH1) 

IF (XND(KC)-1.0)95,34,354 
XNDCKC)=.9999 

CONTINUE 


FIND AVERAGE VELOCITY 


96 


114 


7314 


732 


SUMV=0. 

DO 114 L=i,NNT 
VLECLIFVI*XNDCL) 
SUMV=SUMV4+VL(L) 
VLAV=SUMV/UN 
PLIFT=(2.3769E-3/2.)*VLAVaw2"3.14167%(5./8.)*RPC II) ew 2/444, 
Pe y= =P T 
VOLA46/5.*S-149159*RP CTI] ) eed 
DO 730 MO=1,4 

GO 106731,7323753,754),M0 
DNS=75. 

GO 10 7S> 

DNS=90. 

GO: 10..735 


95 





7S NNS=#109. 
G0.10-735 
734 DNS=125. 
735 PWI=DNS/1728.*VOL 
- PRINT 736, DNS 
736 FORMAT(/18HPARTICLE DENSITY = F6.1) 
PLIOT= PLIFT-PWwI+FLSC 
PMAS=PWT/52.16 
PLAC=PLTOT/PMAS 
IF (PLAC) 842,843,843 
B42 PLAC=0.0 
845 PRINT 110 
110 FORMAT(/22HAVERAGE LOCAL VELOCITY,3X,16HLIFT ON PARTICLE,3SX, 
1 26HNET LIFT FORCE ON PARTICLE,4X,21HPARTICLE ACCELERATION) 
PRINT 120 
120 FORMAT(7X,6H(FT/SEC),7X»22HDUE TO WALL CONSTRAINT, SOX,21HIN VERTIC 
1AL DIRECIIONS) 
PRINT 840,VLAV,PLIFT,PLTOT,PLAC 
840 FORMAT(F20.10-E20.10,5X,F20.10,5X,F20.10/) 
730 CONTINUE 
GO TO 115 
666 PRINT 667 
667 FORMAT(/SHERROR) 
GO TO 668 
115 CONTINUE 
668 RETURN 
END 
SUBROUTINE ROOT(NNT,»MM,PCOZ,YND,XND,NCC) 
DIMENSION YND(70)+XND(70),PCOZ(30),DVX(30),FX(30) 
DO 210 L=i,NNT 
IF(NCC-24)199,199,201 
199 IF CYND(L)-.425)200,201,201 
201 XND(L)=0.1 | 
106 SFX=0. 
SDVX=0. 
DOLOON=1,M™M 
MXL=EMM-N+1 
MX2=MM-N 
EXL=MX4 
FX(N)=PCOZ(N) *XND(L) **MX1 
IF (N-MM)1012102,101 
1014 DVX(N)=EX1*PCOZ(N) *XND(L) #&*MX2 
GO To 1038 
102 DVX(N)=PCOZ(N) 
103 SFX=SFX+FX(N) 
100 SDVX=SDVX+DVX(N) 
RMX=(SFX-YND(L))/SDVX 
IF CABSFCRMX)-.01 )108,108,104 
104 XND(L)=XwD(L)-RMX 
MLM=1000/MM 
IF CXND(L)-2.*"*MLM)106,2049,204 
GO TO 210 


96 





DO GrateCNCCc 471055209, 410 

109 IP(NNE=2)210.210,110 

iq elre COND CI) —1.0)2105 2065206 
205 IF CYND(L)-2.0)204,204,205 
202 XND(L)=.9999 


GO TO 210 


204 3F(L-1)304, 5094, $09 


$04 XND(L)=0.0 
GO TO 210 


S05 XxnD(L)=XNN(L-1) 


GO TO 210 
200) GO. | or <ci10 
210 CONTINUE 

RETURN 

END 

END 


DRAG FORCES: 
Input. 


Program 
program. 


is generated in the same manner as the lift 


Intermediate Results and Output. 


DELT 
XDEL 
poo 
PME 
PWT 
PDAG 


Boundary layer thickness 
at top of boundary layer 
Drag due to pressure force 
Drag Gue ELoudrici ton 
Particle weight 
Horizontal acceleratiog of particle 


SUBROUTINE DRAG(PCOZ,MM,NCC,VJ,RG,ZG,DN» YHVX,VMX) 


DIMENSION 


DP(10),UN(40),RP(10),YGU(40),YGB(490),YL(70),YG(70), 


1 ARC7O) 5 YNDS 70 eX NU7 DO) s VIC 70) PEOZ (S00. Xl 40) X50. 
2 XDEL(5),YMAX(5),PCOX(10),ARF(70) 


DP(1)=.06 
DP(2)=.04 


Hf 





6114 
612 
400 


401 


403 


503 


302 
661 
649 
61 
19 
60 
18 
20 
22 
644 
644 
645 


647 


OP iso = oe 

DP(49)=.91 

DP C5) =20.0 7 

be GO y= te> 

DP(7)=.003 

NN=S4 

RGN=RG/DN 

ZGN=ZG/DN 

PRINT 611, LN 

FORMAT(C22HNOZZLE DIAMETER (FT) = F5.2) 
PRINT 612, VJ 

FORMAT (C26HNOZZLE VELOCITY (FT/SEC) = F7.1i ) 
PRINT 4990,RKGN 

FORMAT(26HNONDIMEN RADIAL DISTANCE = F20.10) 
PRINT 401, ZGN 

FORMAT(28HNONDIMEN HEIGHTH OF NOZZLE = F2Q.10) 
GO.495 | l=197 

PRINT 404, DP(IT) 

FORMAT( 23HPARTICLE DIAMETER-IN. = F2Q-10 ) 
RPCII)V=DPCII)/2, 

UNIT =NN 

TKE=RPCII)/(UNI-“,5) 

YibC= G4 

SUMAR=0,. 

NAT=NN-1 

DO 38 N=L,NAT 

YL(N)=YLC+TK/2, 

AM=N 

XLCN)SSQRTFECRPC II) **2-YL(N) *¥2) 
AR(N)=2.*XLON)¥TK 

ARF (N)=4.*XLO(N) *TK 

IF (N-1)8,5,58 
SUMAR=SUMAR+3.1416/2. *¥ ARF (N) 

EO TG 6 

SUMAR=SUMAR+3.1416"ARF(N) 
YGUCN)=RPCII)+*YLC 

YGBCN) =RPCII)-YLC 

YLC=AM*TK 

GNU=,.00000038745/7.0023769 

Pe CNCC-25)502,505,502 

VI=VMX 

GO TO 700 

IF (ZG/DN-.5)640,640,661 

KZG=ZG/DN 

GO T0(640,641,643,643),K2G 

IF (RGN-.35)60,60,61 

IF (RGN-.75)18.18,19 

IF (RGN-1.2)20,20,22 

VIZFVU*SQRTF CRGNw*2+,.0175*RGN) 

GO TO 700 

VIZEVU*SORTF(-3,2eRGNeH2+5 LIeRGN-1,26) 
GO 10 708 
VIZVU¥SORTE (-2, 68% KGN ee 2+5,48e¥RGN-1. 82) 
GO 10 700 
VIZVU*SORTF(.256eRGNew 2-1, 48 eRGN+2,511) 
GO TO 700 

IF (RGN-,35)69,60,644 


‘TF (CRGN-1.0)645,6495,647 


VIZVI*SORTE (-.0985*RGNee2+4,135"RGN-.241) 
GO 02706 
VIZVU*SQRTFE(.LO77#®RGNw wl, 641 eRGN+1, 5549) 


GO TO 700 
98 





643 
654 
655 


657 


62 
700 


24 
26 
ay. 


28 
977 
900 


99 
S01 
978 


29 
672 


804 
800 
811 
801 
802 
806 
807 
808 
809 
805 
75 


33 
72 


73 
670 
34 
$2 


74 
79 


76 


IF (RGN~.40)62,62,654 

IF (RGN-1.2)655,655,657 
VIZEVU*SQRTF (-. 5 *RGN#© e241, 2"RGN-.S) 
GO TO 709 
VIZVU*SORTFE (0904 #RGNew2- , S1BeRGN+, 744) 
GOTO 700 
VI=SVU*SORTF (625 eRGNe wee, 025eRGN) 
NNT=2*NN-3 

DOZYN=1,NNI 

IF (N-NN) 24, 26,26 

YG(N) ®YGBCNN=N) 

GO. 10-27 

YG(N) =YGUCN-NN4#2) 

IF (N“NN)27,27,28 

ARCN) SARCNNAN) 

GO TO 977 

AR(N) =AR(N-NN#2) 

IF (NCC-26)99,900,900 

AAHSV1I/RG 

YNDCN) =YG(N) *SQRTF (2. *AAA/GNU) /12. 
GO TO 29 

IF (NCC-25)978,301,900 

YND(N) =YG(N) /YHVX 

GO TO 29 

YND(CNI =CYG(N)/(12,*DN))*SQRTF(VU*®DN/GNU) 
CONTINUE 

DO 805 NW=iL,NNT 

XND(NW) 20.0 

IF (NCC-25)808,804,808 

IF CYND(NW)-.15)800,800,801 
MM=8 

DO 811 I2Z=i.5MM 
PCOX(IZ)=PCOZ2(1Z) 

GO TO 808 

MM=4 

DO 80e¢ J=1,MM 
PCOX(I)=PCOZ(1+8) 

DO 807 Kw=i1.MM 

EX1=MM-KW 

FX (KW) =PCOX (KW) *YND(NW) *#EX1 
XNDCNW)=XNDCNW) #FXCKW) 

GO TO 805 

DO 809 KY=1,MM 

EXLEMM-KY+i 
FX(KY)=PCOZ( KY) *®YND(NW) *wEX1 
XNDONW)=XNDCNW) #PXCKY) 
CONTINUE 

DO S2 KC=1,. NN) 

IF (XND(KC) 233,670,670 

IF (KC-1)72,72,73 

VND CKCJ=0;.0 

GO 10-670 

XND( KC) =XNDCKC-1) 

IF (XND(KC)-1,.0)32,52,54 
XND(KC)2.9999 

CONTINUE 

IF (11-1)74,74,80 
DELT=4.6*12./SQRTF(2.*AAA/GNU) 
GO TO 80 

XDEL=0.99 

CALL ROOT(1.MM,PCOX,XDEL,YMAX,NCC) 
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DELT=YMAX(1)*YHVX 
GOo 10. 20 
ZO PEC a4. 99 
CALL RUO1(1,.MM,PCOZ,XDEL,YMAX,WNCC) 
DELT=YMAX(1)*DN/SGRIF (VJI*DN/GNU) #12. 
80.60 TO $2 
$4 DOSSK=1>NNi 
35 VLOCK)=VI*¥XNDCK) 
62 NMID=S(NNI41)/2 
RN=VLONMID) *DPCII)/(12. *GNU) 
LP CRN=100.)74505760, 784 
73% VPURN=100003)752.752,755 
730 CD=S$0./(RN#**, 699) 
GO TO 734 
732 CD=1.494/(RN«¥xw 9396) 
Goa TO. 734 
755 .C0=1.138 
754) Fa 04 
DOS6K=1,NNT 
FD=2.S$769E-3*CDeVL(K) we2wAR(K)/288. 
36 FDI=FDT+FD 
ARS=4-%3.14159*"RP( 11) ««2 
ARD=ARS-SUMAR 
ARD=ARD/ARS 
ACOR=1.+ARD 
CF= .664/SQRTF CRN) 
FFDT=0. 
DO4OM=1,NNI 
FFD=2.$769E-S*CFeVL(M) ew2e ARF (M)/288. 
40 FFDT=FFDI+FFD 
CFB=4.%*CF 
TOL=ABSF(CEB«.00001) - 
CALL CUBE(CFB,TOL,CF83) 
CBD=0.1/CFBS 
FBDT=0., 
DO 42 MD=1,NNT 
FBD=2.3769E-3*CBD*eVL(MD) *w2*AR(MD)/288, 
42 FRDT=FBDT+FRD : 
PRINT $8 
38 FORMAT(/4X,16HPRESSURE DRAG-LB,4X,16HFRICTION DRAG-LB,8X,12HBASE [ 
LRAG“-LB,2X%,24HBOUNDARY LAYER THICKNESS,2X,351HSPHERE APPROXIMATION [ 
2IFFERENCE//) 
PRINT 41, FUT,FFDT,FBDT,DELT,ARD 
41 FORMAT( 3&20.10,2F20.10) 
VOL=4~./35.*5.14159"RP( IL] ) wud 
GS MG(72is 722,725,756) 4M0 
721 UNS=75. 
GO 10, 724 
722 DNS=90. 
GO TO 724 
725° GNS=105. 
GO TO 724 
756 DNS=125. 
724 PWI=DNS/1728.*VOL 
PRINT 725,DNS 
72> FORMAT C/ICAPART TOME SDENSITY¥s= F 6¢%)) 
TOIL=(FDI*FFDT-,707*PWT)*ACOR 
PRINT 620, TOTL 
620 FORMAT(/S54HNET FORCE IN DIRECTION OF MOTION = €E20,10) 
PMAS=PWT/32.16 
PDAC=TOTL/PMAS 


LOO 





IF (PDAC) 842,843,843 
842 PDAC=0.0 
845 PRINT 621, PDAC 
ee teen ACCELERATION IN HORIZONTAL DIRECTION (FT/SEC) 
= Oi) 
720 CONTINUE 
415 CONTINUE 
IF (NCC-24) 668,668,599 
599 CONTINUE 
668 REIURN 
END 


101 





+A 








e) r - 


SPP ESO Sh at 2478s 
h. 023408 « 

SACS a fs TeIAG 4 
JOT THASHHEN)  AMNOFT 
fee. oS" s 
“Unaltaod 

~ Meet 4 Woy 

66 «Bdd(hS*09N 15), 
1wo9 

ARUTLSIR } 

Urs 


7 
r) 




















